classdef NMRGyroscope < handle
    %NMRGYROSCOPE Summary of this class goes here
    %   Detailed explanation goes here
    
    properties
        session
        
        ctr
        rpm
        nmr
        
        plt
        plt_data
        data
    end
    
    methods
        function obj = NMRGyroscope(session)
            %NMRGYROSCOPE Construct an instance of this class
            %   Detailed explanation goes here
            L = log4m.getLogger();
            
            obj.session = session;
            obj.ctr = NMRGyroscope.DeviceController.MainDeviceController(session);            
            obj.rpm = NMRGyroscope.ParametricMagnetometer.RubidiumParametricMagnetometer(session, obj.ctr);
            obj.nmr = NMRGyroscope.NMRSystem.XeNMRSystem(session, obj.ctr);
            
            L.trace('NMRGyroscope', 'NMRGyroscope object created.');
        end
        
        function initialization(obj)
            obj.ctr.SwitchMachineState('init');
            
            obj.rpm.setMainField(); % set main field, use default value
            obj.rpm.GET_WorkingPoint();
            obj.rpm.lockBy();
            obj.rpm.lockBz();
            obj.nmr.DO_FindNMRResonance();
            obj.nmr.GET_OptimalPIDParameters();
            obj.nmr.lockXe129();
            obj.nmr.lockXe131();
        end
    end    
    
    %% Deprecated
    methods
        function obj=getNMRSignal(obj, varargin)
            p=inputParser;
            p.addParameter('PID1_P', -7e-3, @isnumeric);
            p.addParameter('PID1_I', -300e-6, @isnumeric);
            p.addParameter('PID1_D', 0.0, @isnumeric);
            p.addParameter('PID2_P', -1e-3, @isnumeric);
            p.addParameter('PID2_I', -200e-6, @isnumeric);
            p.addParameter('PID2_D', 0.0, @isnumeric);
            p.addParameter('pidFilterBW1', 5.0, @isnumeric);
            p.addParameter('pidFilterBW2', 5.0, @isnumeric);
            p.addParameter('window_time', 13.0, @isnumeric);
            p.addParameter('total_duration', 12.0, @isnumeric);
            p.parse(varargin{:});
            
            ziDAQ('setDouble', '/dev3560/pids/0/p', p.Results.PID1_P);
            ziDAQ('setDouble', '/dev3560/pids/0/i', p.Results.PID1_I);
            ziDAQ('setDouble', '/dev3560/pids/0/d', p.Results.PID1_D);
            ziDAQ('setDouble', '/dev3560/pids/1/p', p.Results.PID2_P);
            ziDAQ('setDouble', '/dev3560/pids/1/i', p.Results.PID2_I);
            ziDAQ('setDouble', '/dev3560/pids/1/d', p.Results.PID2_D);
            ziDAQ('setInt', '/dev3560/pids/0/enable', 1);
            ziDAQ('setInt', '/dev3560/pids/1/enable', 1);
            
            obj.ziSlave.setDemodFilterBW('2', p.Results.pidFilterBW1);
            obj.ziSlave.setDemodFilterBW('3', p.Results.pidFilterBW2);
            
            obj.runPlotter('windowTime', p.Results.window_time, 'totalDuration', p.Results.total_duration);
            obj.data.time=obj.plt_data{1}.time;
            obj.data.frequency129=obj.plt_data{1}.value;
            obj.data.frequency131=obj.plt_data{2}.value;
            obj.data.angleRate=obj.getAngleRate;
                                  
            obj.data.fid.time=obj.data.time;
            obj.data.fid.value=obj.data.angleRate;
            obj.data.fft=obj.fftSignal(obj.data.fid);
                        
            obj.plotData;
            obj.nmr.saveData('NMRSignal');
            
        end
        
        function obj=runPlotter(obj, varargin)
            p=inputParser;
            p.addParameter('windowTime', 13.0, @isnumeric);
            p.addParameter('totalDuration', 12.0, @isnumeric);
            p.parse(varargin{:});
            
            obj.plt=obj.ziSlave.plotterModule;
            obj.plt.reset();
            obj.plt.subscribeChannel('0', {'frequency'});
            obj.plt.subscribeChannel('1', {'frequency'});
            pltData = obj.plt.run('window_time', p.Results.windowTime, 'total_duration', p.Results.totalDuration);
            obj.plt.unsubscribeChannel();
            obj.plt_data{1} = pltData.getNodeData('0', 'frequency');
            obj.plt_data{2} = pltData.getNodeData('1', 'frequency');        
        end
             
                     
        function angleRate=getAngleRate(obj)                     
            angleRate=zeros(1, length(obj.data.time));
            for k=1:length(obj.data.time)
                angleRate(k)=(abs(obj.gamma129)*obj.data.frequency131(k)-obj.gamma131*obj.data.frequency129(k))/(obj.gamma131-abs(obj.gamma129));
            end
        end
        
        function fftData = fftSignal(obj, fidData)
            samplingRate=1/(fidData.time(2)-fidData.time(1)); 
            [fftData.amplitude, fftData.frequency]=periodogram(fidData.value, rectwin(length(fidData.value)), length(fidData.value), samplingRate, 'psd');
            obj.plotFFTData(fftData);
        end
        
        function obj=plotData(obj)
            figure('Name', 'Plotter data');
            subplot(3,1,1)
            plot(obj.plt_data{1}.time, (obj.plt_data{1}.value-mean(obj.plt_data{1}.value))*1e3, obj.plt_data{2}.time, abs(obj.gamma129)/obj.gamma131*(obj.plt_data{2}.value-mean(obj.plt_data{2}.value))*1e3);
            xlabel('Time (s)');ylabel('Frequency shift (mHz)');
            xlim([min(obj.plt_data{1}.time), max(obj.plt_data{1}.time)]);
            legend(sprintf('Xe129, std=%0.4f mHz', std((obj.plt_data{1}.value-mean(obj.plt_data{1}.value))*1e3)),...
                   sprintf('Xe131, std=%0.4f mHz', std((obj.plt_data{2}.value-mean(obj.plt_data{2}.value))*1e3)));
            
            subplot(3,1,2)
            plot(obj.plt_data{1}.time, obj.data.angleRate*1e3)
            xlabel('Time (s)');ylabel('Frequency (mHz)');
            xlim([min(obj.plt_data{1}.time), max(obj.plt_data{1}.time)]);
            legend(sprintf('Angular rate, std=%0.4f mHz', std(obj.data.angleRate*1e3)));
            
            subplot(3,1,3)
            fftData.frequency=obj.data.fft.frequency;fftData.amplitude=sqrt(obj.data.fft.amplitude)*21600; % degree/sqrt(h)
            obj.plotFFTData(fftData, 'legName', 'FFT of angular rate');
        end
        
        function obj = plotFFTData(obj, fftData, varargin)
            p=inputParser;
            p.addRequired('fftData', @isstruct);
            p.addParameter('legName', 'FFT', @ischar);
            p.parse(fftData, varargin{:});
            
%             figure
            loglog(fftData.frequency, fftData.amplitude)
            xlabel('Frequency (Hz)');ylabel('Spectral density (Degree/sqrt(h)');
            xlim([min(fftData.frequency), max(fftData.frequency)]);
            legend(p.Results.legName);
        end
        
    end
end

